home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / SATELITE / ORBIT23.ZIP / ORBITR.C < prev    next >
C/C++ Source or Header  |  1987-01-28  |  8KB  |  293 lines
  1. /* N3EMO Orbit Simulator routines */
  2.  
  3. /* Copyright (C) 1986 Robert W. Berger
  4.    May be used for non-commercial purposes without prior written permission. */
  5.  
  6. #include <stdio.h>
  7. #include <math.h>
  8.  
  9. #define PI 3.14159265
  10. #define PI2 (PI*2)
  11. #define MinutesPerDay (24*60.0)
  12. #define SecondsPerDay (60*MinutesPerDay)
  13. #define HalfSecond (0.5/SecondsPerDay)
  14. #define EarthRadius 6378.16             /* Kilometers */
  15. #define EarthFlat (1/298.25)            /* Earth Flattening Coeff. */
  16. #define SiderealSolar 1.0027379093
  17. #define DegreesPerRadian (180/PI)
  18. #define RadiansPerDegree (PI/180)
  19. #define ABS(x) ((x) < 0 ? (-(x)) : (x))
  20. #define SQR(x) ((x)*(x))
  21.  
  22. #define RefYear 79
  23. #define RefSidereal 0.27518504
  24. #define RefDay 0                /* Day of week */
  25.  
  26. char *MonthNames[] = { "Jan","Feb","Mar","Apr","May","Jun","Jul",
  27.                         "Aug","Sep","Oct","Nov","Dec" };
  28.  
  29. int MonthDays[] = {31,28,31,30,31,30,31,31,30,31,30,31};
  30.  
  31. char *DayNames[] = { "Sunday","Monday","Tuesday","Wednesday","Thursday",
  32.                         "Friday","Saturday"};
  33.  
  34. /* Solve Kepler's equation                                      */
  35. /* Inputs:                                                      */
  36. /*      MeanAnomaly     Time since last perigee, in radians.    */
  37. /*                      PI2 = one complete orbit.               */
  38. /*      Eccentricity    Eccentricity of orbit's ellipse.        */
  39. /* Output:                                                      */
  40. /*      TrueAnomaly     Angle between perigee, geocenter, and   */
  41. /*                      current position.                       */
  42.  
  43. double Kepler(MeanAnomaly,Eccentricity)
  44. register double MeanAnomaly,Eccentricity;
  45.  
  46. {
  47. register double E;              /* Eccentric Anomaly                    */
  48. register double Error;
  49. register double TrueAnomaly;
  50.  
  51.     E = MeanAnomaly;    /* Initial guess */
  52.     do
  53.         {
  54.         Error = (E - Eccentricity*sin(E) - MeanAnomaly)
  55.                 / (1 - Eccentricity*cos(E));
  56.         E -= Error;
  57.         }
  58.    while (ABS(Error) >= 0.000001);
  59.  
  60.     if (ABS(E-PI) < 0.000001)
  61.         TrueAnomaly = PI;
  62.       else
  63.         TrueAnomaly = 2*atan(sqrt((1+Eccentricity)/(1-Eccentricity))
  64.                                 *tan(E/2));
  65.     if (TrueAnomaly < 0)
  66.         TrueAnomaly += PI2;
  67.  
  68.     return TrueAnomaly;
  69. }
  70.  
  71. GetSubSatPoint(EpochTime,EpochRAAN,EpochArgPerigee,Inclination,
  72.         Eccentricity,RAANPrecession,PerigeePrecession,Time,TrueAnomaly,
  73.                 Latitude,Longitude)
  74.  
  75. double EpochTime,EpochRAAN,EpochArgPerigee,Inclination,Eccentricity;
  76. double RAANPrecession,PerigeePrecession,Time;
  77. double TrueAnomaly,*Longitude,*Latitude;
  78. {
  79.     double RAAN,ArgPerigee;
  80.     int i;
  81.  
  82.     ArgPerigee = EpochArgPerigee + (Time-EpochTime)*PerigeePrecession;
  83.  
  84.     *Latitude = asin(sin(Inclination)*sin(ArgPerigee+TrueAnomaly));
  85.  
  86.     RAAN = EpochRAAN - (Time-EpochTime)*RAANPrecession;
  87.  
  88.     *Longitude = -acos(cos(TrueAnomaly+ArgPerigee)/cos(*Latitude));
  89.     if ((Inclination > PI/2 && *Latitude < 0)
  90.           || (Inclination < PI/2 && *Latitude > 0))
  91.         *Longitude = -*Longitude;
  92.     *Longitude += RAAN;
  93.  
  94.     /* Convert from celestial to terrestrial coordinates */
  95.     *Longitude -= PI2*(Time*SiderealSolar + RefSidereal);
  96.  
  97.     /* Make West be positive */
  98.     *Longitude = -*Longitude;
  99.  
  100.     /* i = floor(Longitude/2*pi)        */
  101.     i = *Longitude/PI2;
  102.     if(i < 0)
  103.         i--;
  104.  
  105.     *Longitude -= i*PI2;
  106. }
  107.  
  108.  
  109. GetPrecession(SemiMajorAxis,Eccentricity,Inclination,
  110.         RAANPrecession,PerigeePrecession)
  111. double SemiMajorAxis,Eccentricity,Inclination;
  112. double *RAANPrecession,*PerigeePrecession;
  113. {
  114.   *RAANPrecession = 9.95*pow(EarthRadius/SemiMajorAxis,3.5) * cos(Inclination)
  115.                  / SQR(1-SQR(Eccentricity)) * RadiansPerDegree;
  116.  
  117.   *PerigeePrecession = 4.97*pow(EarthRadius/SemiMajorAxis,3.5)
  118.          * (5*SQR(cos(Inclination))-1)
  119.                  / SQR(1-SQR(Eccentricity)) * RadiansPerDegree;
  120. }
  121.  
  122. GetBearings(SiteLat,SiteLong,SiteAltitude,
  123.         SatLat,SatLong,Radius,Azimuth,Elevation,Range)
  124. double SiteLat,SiteLong,SiteAltitude,SatLat,SatLong,Radius;
  125. double *Azimuth,*Elevation,*Range;
  126. {
  127.     double SinSiteLat,sinSiteLong,SinSatLat,SinSatLong;
  128.     double CosSiteLat,cosSiteLong,CosSatLat,CosSatLong;
  129.     double CosBeta,SinBeta;
  130.     double LongDiff;
  131.  
  132.     SinSiteLat = sin(SiteLat); sinSiteLong = sin(SiteLong);
  133.     CosSiteLat = cos(SiteLat); cosSiteLong = cos(SiteLong);
  134.     SinSatLat = sin(SatLat); SinSatLong = sin(SatLong);
  135.     CosSatLat = cos(SatLat); CosSatLong = cos(SatLong);
  136.  
  137.     LongDiff = SiteLong - SatLong;
  138.     CosBeta = SinSiteLat*SinSatLat+CosSiteLat*CosSatLat*cos(LongDiff);
  139.     SinBeta = sqrt(1-SQR(CosBeta));
  140.  
  141.     *Azimuth = acos((SinSatLat- SinSiteLat*CosBeta)/CosSiteLat/SinBeta);
  142.  
  143.     if (LongDiff < -PI)
  144.         LongDiff += PI2;
  145.     if (LongDiff > PI)
  146.         LongDiff -= PI2;
  147.  
  148.     if (LongDiff < 0)
  149.         *Azimuth = PI2 - *Azimuth;
  150.  
  151.     /* Convert to geocentric radius */
  152.     SiteAltitude += EarthRadius*(1-EarthFlat/2+EarthFlat/2*cos(2*SiteLat));
  153.  
  154.     *Elevation = atan((Radius*CosBeta-(SiteAltitude))
  155.                         /(Radius*SinBeta));
  156.  
  157.     *Range = sqrt(SQR(Radius) + SQR(SiteAltitude)
  158.                     -2*Radius*SiteAltitude*CosBeta);
  159. }
  160.  
  161.  
  162. SPrintDate(Str,Day)
  163. char *Str;
  164. {
  165.     int Month,Year,DayOfWeek;
  166.  
  167.     DayOfWeek = (Day-RefDay) % 7;
  168.  
  169.     Year = RefYear;
  170.  
  171.     while (Day <= 0)
  172.         {
  173.         Year -= 1;
  174.         if (Year%4 == 0)
  175.             Day += 366;
  176.          else
  177.             Day += 365;
  178.         }
  179.  
  180.     while ((Year%4 == 0 && Day > 366) || (Year%4 != 0 && Day > 365))
  181.         {
  182.         if (Year%4 == 0)
  183.             Day -= 366;
  184.          else
  185.             Day -= 365;
  186.         Year += 1;
  187.         }
  188.  
  189.     if (Year%4 == 0)
  190.         MonthDays[1] += 1;      /* Leap year */
  191.  
  192.     Month = 0;
  193.     while (Day > MonthDays[Month])
  194.         {
  195.         Day -= MonthDays[Month];
  196.         Month += 1;
  197.         }
  198.  
  199.     sprintf(Str,"%s  %d %s %d",DayNames[DayOfWeek],Day,
  200.                 MonthNames[Month],Year);
  201.  
  202.     if (Year%4 == 0)
  203.         MonthDays[1] -= 1;      /* Leap year */
  204. }
  205.  
  206.  
  207. SPrintTime(Str,Time)
  208. char *Str;
  209. double Time;
  210. {
  211.     int day,hours,minutes,seconds;
  212.  
  213.     day = Time;
  214.     Time -= day;
  215.     if (Time < 0)
  216.         Time += 1.0;   /* Correct for truncation problems with negatives */
  217.  
  218.     hours = Time*24;
  219.     Time -=  hours/24.0;
  220.  
  221.     minutes = Time*MinutesPerDay;
  222.     Time -= minutes/MinutesPerDay;
  223.  
  224.     seconds = Time*SecondsPerDay;
  225.     seconds -= seconds/SecondsPerDay;
  226.  
  227.     sprintf(Str,"%02d%02d:%02d",hours,minutes,seconds);
  228. }
  229.  
  230. PrintDate(OutFile,Day)
  231. FILE *OutFile;
  232. {
  233.     char str[100];
  234.  
  235.     SPrintDate(str,Day);
  236.     fprintf(OutFile,"%s",str);
  237. }
  238.  
  239. PrintTime(OutFile,Time)
  240. FILE *OutFile;
  241. double Time;
  242. {
  243.     char str[100];
  244.  
  245.     SPrintTime(str,Time);
  246.     fprintf(OutFile,"%s",str);
  247. }
  248.  
  249.  
  250. /* Get the Day Number for a given date. January 1 of the reference year
  251.    is day 1. Note that the Day Number may be negative, if the sidereal
  252.    reference is in the future.                                          */
  253.  
  254. double GetDay(Year,Month,Day)
  255. double Day;
  256. {
  257.     int TmpYear;
  258.  
  259.     TmpYear = Year;
  260.  
  261.     while (TmpYear > RefYear)
  262.         {
  263.         TmpYear -= 1;
  264.         if (TmpYear%4 == 0)
  265.             Day += 366;
  266.          else
  267.             Day += 365;
  268.         }
  269.  
  270.     while (TmpYear < RefYear)
  271.         {
  272.         if (TmpYear%4 == 0)
  273.             Day -= 366;
  274.          else
  275.             Day -= 365;
  276.         TmpYear += 1;
  277.         }
  278.  
  279.     if (Year%4 == 0)
  280.         MonthDays[1] += 1;      /* Leap year */
  281.  
  282.      while (Month > 1)
  283.         {
  284.         Month -= 1;
  285.         Day += MonthDays[Month-1];
  286.         }
  287.  
  288.     if (Year%4 == 0)
  289.         MonthDays[1] -= 1;      /* Leap year */
  290.  
  291.     return Day;
  292. }
  293.